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Abstract 

Problems in R 3 are addressed where the scalar potential of an associated vector field satisfies Laplace's 
equation in some unbounded external region and is to be approximated by unknown (point) sources con- 
tained in the complimentary subregion. Two specific field geometries are considered: R 3 half-space and 
the exterior of an R 3 sphere, which are the two standard settings for geophysical and geoexploration 
gravitational problems. For these geometries it is shown that a new type of kernel space exists, which is 
labeled a Dirichlet-integral dual-access collocation-kernel space (DIDACKS) and that is well suited for 
many applications. The DIDACKS examples studied are related to reproducing kernel Hilbert spaces and 
they have a replicating kernel (as opposed to a reproducing kernel) that has the ubiquitous form of the in- 
verse of the distance between a field point and a corresponding source point. Underpinning this approach 
are three basic mathematical relationships of general interest. Two of these relationships — corresponding 
to the two geometries — yield exact closed-form inner products and thus exact linear equation sets for 
the corresponding point source strengths of various types (i.e., point mass, point dipole and/or point 
quadrupole sets) at specified source locations. The given field is reconstructed not only in a point 
collocation sense, but also in a (weighted) field-energy error-minimization sense. 
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kernels, fundamental solutions, point sources, multipole, potential theory 

AMS subject classifications. 35J05, 31B10, 86A22, 65D05 

1 Introduction 

The goal of this article is to set forth a mathematical framework for the approximation of R 3 harmonic 
fields in unbounded domains by point sources contained inside the complimentary region. The proposed 
Dirichlet-integral dual-access collocation-kernel (DIDACK) approach has a mathematically and physically 
well-motivated underpinning. The associated space (DIDACKS) has certain similarities to reproducing 
kernel Hilbert space (RKHS), but is distinct from it. Two concrete R 3 geometries are considered: (A) The 
harmonic field region consisting of half-space (denoted fix). (B) The harmonic field region consisting of the 
exterior of a sphere (denoted Qq). Within this geometric context, the developed formalism easily handles 
various combinations of diverse types of point sources (such as point masses, point mass dipoles or point 
mass quadrupoles) ; moreover, for a set of specified source locations the formalism yields closed- form linear 
equation sets that simultaneously minimize the volume integrals of (weighted) field energy densities [i.e., 
(weighted) Dirichlet integrals]. 
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Techniques introduced here can be either applied directly or adapted for use in many mathematical and 
physical areas. Examples on the mathematical side include potential theory, point source elliptic bound- 
ary value modeling (i.e., method of fundamental solutions), fast multipole techniques, radial-basis function 
techniques, RKHS techniques, geophysical collocation (GC) techniques and standard energy minimization 
based techniques (such as Galerkin and Raleigh- Ritz based approaches). Also see [21J . Examples on the 
physical side include geoexploration (where gravity is often used to locate oil or other minerals), geophysics, 
magnetostatics (for a survey of these three areas see [TTJ [24] [22J, [4] ) , general electromagnetic source analysis, 
electrostatics and hydrostatics (for the physical significance of these three areas, and of field energy with 
regards to them, see [3]). Also somewhat similar mathematical structures arise in many other areas such as 
biophysical or biomedical engineering {where, for example, electric dipole models are used in electrocardio- 
graphical (ECG) modeling and current dipole models are used in electroencephalography (EEG) [5]}. 

First, given the diversity just noted, it is clearly desirable to not confine potential readership to one 
specialty or another; hence, the required mathematical background for the core part of the article is limited 
to classical R 3 potential theory [TU [TU] and a basic understanding of functional analysis, as utilized in 
approximation theory. This core part is primarily mathematical in content and makes up the first four 
sections. Another, perhaps stronger, motivation exists for limiting the level of this core material — there is 
often a given mathematical level that certain research results tend to emerge from and the same results may 
not be apparent at either a higher or lower level of mathematical sophistication. The first four sections, as 
well as the geophysical applications setting given in Section^ are aimed at this default level of sophistication. 
While a knowledge of neither RKHS nor GC is a necessity in these first [5] sections, a better understanding of 
both these areas is presupposed in Section [6] which describes the connections of DIDACKS theory to other 
kernel based approaches. GC, which is often called least-squares collocation, is a RKHS based approach that 
differs from standard collocation techniques utilized by applied mathematicians in several relevant ways and, 
as such, may not be familiar to many readers. Moritz [18j provides a readily accessible treatment of both 
RKHS and GC theory; nevertheless, the thrust of Section [5] should be directly accessible to those familiar 
with only standard RKHS theory {[16], [1] or [26] (which contains PQ)}. 

Next consider the applications setting and why it was chosen. While diverse candidate application areas 
exist, the approach adopted here is to select one particular representative problem area, namely geophysical 
gravitation, and discuss it thoroughly. In addition to Laplacian inverse source theory, geophysical gravitation 
has two separate problem categories that point sources can be used in: field modeling and field estimation. 
For gravitational modeling problems the field is assumed to be known throughout the region of interest and 
a more compact, but accurate, representation is desired. For estimation problems it is assumed that values 
are known accurately at a certain number of points in the field region and that one wants to predict gravity 
values over some part of this field region. Besides combinations of point sources, other techniques (such as 
GC) exist for treating gravitational field modeling and field estimation problems. In a wider mathematical 
context, since the kernels involved satisfy a generalized collocation property, DIDACKS modeling and esti- 
mation approaches described here can be considered harmonic interpolation and extrapolation techniques, 
respectively. With minor differences these gravitation (G) problems use the same standard notation and 
techniques employed in electrostatics, where E is the field [TO] , and so this arena should be readily accessible 
to all interested applied mathematicians, physicist and engineers. (Within this electrostatic context the for- 
malism developed here handles multiple types of various point sources — such as point charges, electrostatic 
dipoles or electrostatic quadrupoles.) For consistency the notation used in the first four sections is also 
specialized to the gravitational setting, but in these sections the intrusion of this setting is minimal. 

Geophysical gravitation as an applications setting was chosen for the following five reasons: (a) It has 
a strong overall historical association with potential theory, (b) It has an easily understood notation and a 
readily accessible mature literature, (c) It provides a family of significant and challenging problems (that 
includes representative geophysical inverse source problems), (d) DIDACKS theory was developed to handle 
gravitational problems and all of the author's direct numerical experience with it is in this arena, (e) GC 
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and its extensions are the family of approaches that are mathematically closest to DIDACKS theory. A 
discussion of point (a) is well outside the scope of the present article, but it is worth noting in this context 
that Gauss played a very significant role in the history of geophysics [9, p. 1]. With regards to point (b), as 
well (c) , see US! ■ (To gauge the advance of physical geodesy compare and contrast [9] with [5] , upon which 
it is based.) With regards to (c) and (d), since no previous geophysical background is presupposed in the 
first four sections of the paper, the relevant overall geophysical and point mass modeling context is supplied 
in Section [5l Since GC is the most commonly used approach for regional gravitational data processing and 
estimation it is briefly described at the end of Section O The pertinent differences between GC and standard 
collocation techniques utilized by applied mathematicians arc also described in Section [5j Finally, Section [5] 
considers the mathematical connections implied by (e). 

The DIDACKS approach was conceived and is best understood on its own merits as a self-contained 
mathematical theory that is independent of the geophysical connections just indicated and it is this exposition 
that occupies the first four sections. Underpinning this mathematical side of the DIDACKS approach are 
three intrinsically interesting relationships that will now be briefly surveyed. First consider the notation 
employed. For overall accessibility, for consistency with the geophysical and electrostatic literature and to 
avoid various possible notational conflicts preference is given to a pedestrian but unambiguous notation: 
Cartesian coordinates are used and overset arrows employed to denote K 3 vectors, X = (x, y, z) T £ R 3 (for 
superscript T := transpose), while for n / 3 n-dimensional vectors and matrices are denoted by lower and 
upper case bold letters, respectively. Ro is used to denote the radius of the sphere associated with Q that 
is centered over the origin: Slo : = {A € M 3 | |A| > Ro}- Likewise for the half-space case, the origin is chosen 
in the plane d^i := {X e M 3 | z = 0} so that := {X e K 3 | z > 0}. (Observe that the overall shape of 
the subscripts here match that of the associated boundary.) Frequently these two settings will be denoted 
by a subscript j (for example, flj), where j = or 1 is always implied. 

Temporarily leaving aside the issue of admissible functions, these relationships can be compactly stated in 
terms of a Dirichlet integral over some connected but possibly unbounded region f2, which are usually denoted 
by T)[v, w] — fff^ Vu • Vro dV for admissible harmonic functions v(X) and w(X), or in terms of the more 
inclusive concept of a weighted Dirichlet integral for the region f2 denoted by D[u, w, ^, Q], where /i = fi(X) 
is the weighting function so that D{v, w, fi, :— fff n Li^v ■ Vic dV. Clearly D[u, w, 1, ft] — D[v, w]. Let 
l^ 1 :— 1/|A — X'\, where X e Clj and X' is in the corresponding closed source region := f^. C f2'- := 
compliment of j. (By convention, generally primed variables occur in Qj and unprimed ones in fij.) Then 
the first two relations give the replication (or generalized collocation) properties of the DIDACKS kernel 
I- 1 : 

(1) D[w, r\ 1, fit] = 2irw(x', y',-z') 

and 

(HO dk r\ mo, no] = 2^ \P\ w(P)/r 2 q 

with fiQ = 1/r ( r := \X\) and where 




(Eb) P = 

Finally the third relationship ties the unweighted Dirichlet integral over f2 to the weighted integral given 
on the left hand side (LHS) of and can be written as 

@b) D[«, w, 1, fi ] = Ro • D[«, w, mo, n ] + (2ttR ) • {v, «j) , 
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where the surface inner product on the right hand side (RHS) here is defined as 



(0b) («, w) a := (1/4tt) / / v(r, 9, 0) «,(r, 0, 0) da 

and where, as in [8j and [9], a and da have the following meaning when associated with the integral of f(X) 



2tt 7T 

(4! // f(r,9, ( f>)da := J J [f(r, 6, cj>)\ 

d>=0 6=0 



sin 9 d6d<f) 

r=R 



for standard spherical coordinates r, 9, <fi. [Occasionally the limits implicit in ([4]) will be stated explicitly for 
emphasis.] Due to the way r dependence enters in ((4]) it can be used to derive expressions that are otherwise 
not obvious, as will be apparent in the sequel. 

Clearly ([1]) and ([2k) give a means of performing closed-form inner products based on the Dirichlet 
integral, while (J3^i) links the weighted Dirichlet inner product to the unweighted one over fio- The general 
approximation and functional analysis framework for these relationships is given in Section[2j The derivation 
of |T]) is given in Section [3] and that of ([2k) and ((3ja,) in Section [4] In order to understand their connection to 
the applications described in Section [5] and to put them in proper historical context, next briefly consider the 
history of these relationships. Originally ((2k) and ((3k) were discovered in a different form by the author in 
the early 1980's — namely ([2"T|) and (|2"2"|) . where the "integral norm" introduced in Section[2]is used in place of 
the weighted Dirichlet integral. The history of these relationships in this form and of their application from 
inception to the mid-1990's can be found in [20] . The germane part of this internal history is summarized and 
updated in Sectional As an aside, although fl} is not explicitly mentioned in [20j . the derivation in Section[3] 
of it is the one originally found by the author in the early 1980's. The direct relationship of the integral norm 
to the weighted Dirichlet norm (|3"2"]) is new. Also unless otherwise noted, the presentation itself (including 
various terms and concepts) is completely new here. Finally, while the general DIDACKS technique is a 
synthesis of several results that seem to have been overlooked by the broader scientific community, for any 
mathematical approach that directly touches on harmonic analysis, RKHS theory and gravimetric inverse 
source theory either precedents or specialized parallel lines of development would seem to be a necessity. The 
known ones for R 3 DIDACKS theory are addressed in Section [6] In one way or another all of the discussions 
of Section [5] pertain to the second DIDACKS relation in weighted Dirichlet form, ([2k). As discussed there, 
((2^) should clearly be viewed as going back to Krarup in [14) , since he derived it there in a directly equivalent 
form. Aside from the author's work, there are no known instances of ([1) and ([3k), or for that matter for the 
second DIDACKS relation expressed in the integral norm form (j27jl . 

To motivate a further exploration of the paper's scope and limits it is useful to envision the possible 
reactions of four typical classes of readers to the above relations. First, consider an applied mathematician, 
physicist, engineer or other scientist who may be familiar with the material in the first four chapters of [10j . 
but is not yet a seasoned practitioner and may benefit by consulting [16], [18] or [12]. (The physical and 
historical importance of Dirichlet integrals and of the associated Dirichlet principle may be obtained from 
other sources [3] [T71 [6].) This reader may find all three of the above field energy expressions somewhat 
surprising: ([T]) and ([2k) because they allow for the closed form evaluation of volume integrals and ([3k) 
since it allows for an unusual reexpression of the volume field energy. This reader may also observe that 
([1]) and ([2^) appear to have some connection to Green's functions and the method of images [10] and may 
recognize ([2b ) as the coordinate transformation part of a Kelvin transformation |12[ I13j . These particular 
connections arise from the nature of Green's functions for flj , but the most direct explanation requires some 
knowledge of RKHS theory. First, as noted in [3], the existence of a closed- form reproducing kernel occurs 
when closed-form expressions for both Neumann and Dirichlet Green's functions exist. Second, for the cases 
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studied here a dual-action collocation kernel (DACK) arises from the result of a reflection, (JXJ) , or Kelvin 
transformation, IPJd), applied to a reproducing kernel of the right form for the relevant geometry. 

Second, while the possible reactions of any number of specialists might also be examined, consider a 
reader who has a particular interest in integral kernels or RKHS theory. For various reasons this reader may 
also find the above relationships somewhat surprising. This reader might, for example, observe that i^ 1 plays 
the role of a kernel and then recall that a symmetric reproducing kernel (SRK) of the form \X — Y^ 1 for X 
and Y both in the same region cannot exist since a SRK must be bounded and this kernel is not. Here X' 
is a fixed interior point and X is in the exterior region so I is bounded, which is a very different situation 
and implies a change of perspective. This in itself could raise further questions since although DACKs are 
not reproducing kernels they have some properties analogous to them. In the previous paragraph the kernels 
studied here were linked to reproducing kernels, but it is unclear whether minimum norm DACKs can arise 
in other ways. Moreover, consideration of DACKs as a separate class of kernels also raises the question as to 
how they fit into our current overall understanding of kernel structures. These and other issues of a general 
nature are outside the scope of the present article, but some specific topics that might interest the second 
type of reader, such as the exact definition of the replication and generalized collocation property mentioned 
above, are addressed in Section [2] 

Third, consider a reader who is very applications and results oriented. Such a reader may be disappointed 
to discover that there is not a table containing numerical examples and results; however, a discussion of point 
mass and point dipole DIDACKS results and their associated applications settings can be found in Section[5l 
where global non-linear least squares (NLLSQ) results are emphasized. Due to the variety and nature of 
regional gravity data, as well as other issues [20] . no known easily replicated example provides generic 
benchmark results, which is normally an expectation for these tabulated examples. Moreover, point source 
fitting problems are part of a general class of problems that are "notoriously" ill-conditioned and problematic 
[H pp. 214-222] so that each problem encountered should be tackled on its own terms, which means that one 
or two simple table examples cannot serve to provide adequate implementation guidance. Unfortunately, 
a thorough discussion of associated implementation strategies is outside the scope of the present article. 
Also observe that a concrete example provides a replication check that can serve as a consistency test for 
implementors, but this point is largely superfluous here since the DIDACKS approach exhibits the generalized 
collocation property with respect to point sources and thus, when implemented correctly, any point mass 
or dipole fit replicates the point field data that was used to produce it in the first place to within allowed 
round-off error and thus any implementation serves as its own self-consistency test. 

Fourth and finally, consider a reader whose primary interest is in the theory and application of Laplacian 
inverse source theory. Since there are many shared implementation pitfalls common to both point source field 
modeling and inverse source estimation problems, the comments just made in the last paragraph are also 
relevant in this context. While specific mathematical tools and implementation strategies are not discussed 
here, readers with solid applications experience should be able to make direct use of the formalism presented. 
These readers may also be interested in the topic of continuous parameterized distributions, which is raised 
in Section [2] Finally, it is worth noting that other source region shapes can be entertained within the 
contexts of the two considered geometries since the only real requirement is that source regions be bounded 
and contained within the compliment of the unbounded harmonic field region. 

2 Generalized Linear-Least Squares Setting 

This section addresses the generalized least-squares (GLLSQ) plan of approach and the associated functional 
space setting. There are numerous approaches closely aligned to the GLLSQ one adopted here — such as the 
Galerkin and Raleigh-Ritz based techniques mentioned in the introduction; however, the acronym GLLSQ is 
introduced to imply an implicit change of perspective. In particular, connections to generalized collocation, 
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as discussed later in the section and explicitly formalized by the GLLSQ collocation condition, are implied 
as well as an approach that is distinct from the usual LLSQ ones where sampling and discretization are 
introduced. Connections to GC are also implied. 

Both LLSQ and GLLSQ approaches minimize some cost function <f>' = \\v — w\\ 2 , where for the problems 
of interest v(X) is a point source potential model and w(X) is some given canonical (or truth) reference 
potential. For a point mass fit with Nk point masses v has the form 

N k 

(5) v(X) = G^ ""■ 



fi \ X ~ X 'k\ 



where G is the Newtonian gravitational constant w 6.6742 x 10 11 m 3 s 2 Kg 1 [tj p. 3]. As previously noted, 
X £ flj and X' k S VLs j , which is a bounded and closed subregion of the open region f^-, so that the (kernel) 

basis functions occurring in ([5]) are always bounded. Further X' k , ^ X' k for all k! ^ k is always assumed. 

Five conventions are adopted here. First, in physics texts the potential function v is interpreted as 
potential energy and ([3]) has a negative sign since all gravitational bodies attract and the resulting force 
is given by the negative of the gradient of the potential. In geophysics these sign conventions are different 
and consistent with (0, but in either case this should cause little difficulty since in fitting problems all 
that is required is that the overall sign conventions for v and w be consistent. Second, it is assumed that 
gravitational force is always acting on a unit test mass [9l p. 4] and thus it will be treated as having 
the units of acceleration [9l p. 45]. Third, physical geodesists distinguish between gravity field quantities, 
which include the effects of the Earth's rotation, and gravitational quantities like (O, which do not [9l p. 
44]. This is a distinction physicists generally do not make since rotational effects can be easily tracked and 
accounted for as required. The physicist's lead is followed here and this distinction is ignored. Fourth, 
both positive and negative masses will be considered a possibility since this is the usual convention adopted 
in point mass fitting approaches. Specifically, for gravity modeling and estimation problems each mk can 
clearly be viewed as a mathematical parameter that can assume either sign. This convention also allows for 
the ready adaptation of material developed here to other areas where both signs can occur. (Even regional 
geoexploration inverse mass density estimation problems can be handled by assuming that all smoothed 
density estimates are with respect to an average or ambient density.) Fifth, it is useful to introduce scaled 
versions of the above potential functions in order to absorb the factor of G: V = v/G and W = w/G. Thus 
the cost function to be minimized becomes (with $ := <j>'/G 2 ): 

(6) <Z> = \\V-W\f = \\V\f-2(V, W) + \\W\\ 2 . 

(Notice that scaling a cost function leaves the minima unchanged.) 

Next consider the philosophy behind the norm selection process. As discussed later in Section [S] the 
minimization philosophy of matching the observations as closely as possible has generally been chosen for 
point mass fitting problems. This philosophy is, however, not necessarily sound in all or even most cases. 
For modeling problems a reference model, which is assumed to be accurate, is given and one wishes to 
match this reference as closely as possible in some physical sense. Here the desire is to minimize the possible 
error differences that will result when this given reference model is replaced by a new point mass (or point 
source) model, which invariably occurs in some sort of software emulation of a physical situation. Thus 
instead of "matching the observables as closely as possible," a sounder strategy is to "minimize the type 
of errors that will lead to the greatest errors in the end product." From Newton's second law, since these 
end-product errors here are most often the direct result of gravity errors it is clear that the difference in 
the given gravity reference field and the developed point mass gravity model should be minimized, say in a 
squared residual sense at a large number of appropriate sample points. As this distribution of sample points 
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becomes uniformly dense over the entire global region of interest the following key integral condition results: 



(7) Minimize ® = JJJ - ^W\ 2 dV = D[W - V, W - V, 1, fi,-] . 

Temporarily leaving aside fundamental issues, such as how to turn the RHS of ([7|) into a proper norm 
structure, consider the general form of the linear equation sets that result from minimizing this type of cost 
function. For concreteness consider the minimization process in R 3 half-space (Oi). Since the RHS of Q is 
already proportional to the field energy, it is natural to consider the half-space (j = 1) energy norm: 

i rrr,*„ *„,. a 



(8) \\ V -w\f Ei := — 1 1 1 \VV-VW\* dV 
(a factor of 87r has been inserted since it often occurs for various field energy expressions in appropriate 

2^2 

units). Because \\V — W\\ E = \\V\\ E + || W|| E - 2(V, W) E , the energy inner-product 

(y, W) Ei := D[V, W, 1, fii]/87r 

is also needed. In particular if V is specified through |5]), with if. :— \X — Xi\, and if W is an appropriate 
reference field, then 

N k N k N k 

(9) \\V-Wf Ei = \\W\f Ei - 2gm fc (C W) Ei +J2 E m k m k ,(q\r k }) Ei . 

k=l fc=lfc'=l 

Taking the partial of Q with respect to rrifc" (for k" — 1, 2, 3, ... , Nk), setting the result to zero, then 
dividing by two yields a linear equation set that can be easily inverted for the mass values, provided that 
{(■fZ, </>) E can be easily computed for (f> = W and 4> = l/lk- The relationship which makes this possible is 
(JTJ) . By introducing TJ^fc/ = (l^ 1 ,£~ k , 1 ) E and Ak — {W,l~^ 1 ) E the linear equation set can be written as 

N k 

(10) ^2 T k,k' m k > = A k . 

k'=l 

For the spherical exterior matters proceed in much the same fashion except that a weight function, 
Ho = 1 jr, must be introduced into ([7]) . Not only is this weighting required to turn i^ 1 into a replicating kernel, 
but since regions closer to the Earth's surface are normally of greater interest for geophysical applications 
than regions further away it is also desirable. 

Obviously with regards to applications (|10|) is pivotal and, as such, warrants at least an informal ex- 
amination. Before proceeding it is useful to clarify the differences between reproducing, replication and 
generalized collocation kernels. When a kernel, like i^ 1 , is not symmetric since its arguments are in different 
domains and there is a relationship such as (JTJ) or (J^i) that allows for closed-form inner-product expressions 
it will be called a replication kernel and be said to have the point replication property. Since reproducing 
kernels are necessarily symmetric they cannot be considered replicating kernels so this terminology distin- 
guishes replication and reproducing kernels. Alternatively, the term generalized collocation property will be 
taken to be a generalization of a point data and/or collocation matching condition and as such includes not 
only the possibility of reproducing kernels, but also replicating kernels. As discussed below, it also allows 
for the possibility that resulting inner products may be obtainable by numerical means (after assuming an 
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underlying replication property also holds) — so long as the resulting inner products, Ak, occuring on the 
RHS of (fTU| can be matched. (The Ak's may also represent empirically obtained data.) 

As a concrete example of situations where numerical integration frequently enters consider fits based on 
the continuous analog of ([S]) where the potential is due to some parameterized density function p: 

Here a = ipt\, ct-2, 03, otN k ) T ■ If p consists of a linear superposition of density basis functions tpk 
{i.e., p = X)fe=i a kipk{X')} then minimizing \\W — V\\ 2 yields a linear equation set similar to (fTU|) . where 
the Ak's and 2fc fc/'s must generally be computed numerically; however, when |T]) or (J^i) is used then great 
simplifications result and continuous distributions are tractable. Besides parameterized volume distributions, 
parameterized surface and line distributions are also obviously possible. One interesting choice for volume 
density basis functions tpk is the use of finite element method (FEM) basis functions. Thus consider the case 
where qj. are taken to be a set of node points over fis and the ipk(X') are chosen to be a set of localized FEM 
basis functions with the property that V'fc'( < Zfc) = 1 if fc' = & and V>fc'(<7fc) = otherwise. The ak determined 
directly from the analog of (fT0|) then represent the density strengths at the node points q{. . 

Any additional structure that can help to clarify the possibilities inherent in (TlO|) is desirable. Towards 
that end briefly consider a suggestive theorem from RKHS theory When a reproducing kernel, with specified 
kernel points Qk G Q replacing the values of X' k in (|5|), is used for a basis functions expansion of V that 
is analogous to (O, and when © is replaced by the corresponding cost function based on the reproducing 
kernel norm, then minimization of this cost function results in a closed-form linear equation set that is just 
like (fTOj) — except that the inner products corresponding to Ak take on the simpler form W(Qk)- This type 
of reproducing kernel fit also satisfies a minimum norm collocation property: the function with the smallest 
associated norm that matches the prescribed data set [i.e., the values W(Qk)} is the one which results from 
solving the analog of (JTUJ) [18, pp. 207-220] . This minimum norm property is a well known functional analysis 
result |16| and it insures that a reproducing kernel fit will simultaneously match the given point data and 
minimize both and ||V — W\\ for the associated norm. This fact is of interest here since it strongly 
suggests that if a replicating kernel expansion for V is used in (fTt)|). then generally any specified values for 
Ak are recovered and that D[V, V, 1, fli] or D[V, V, po, Ho] also is simultaneously minimized. While this 
may fail to happen due to auxiliary restrictions placed on the basis functions or on the overall space of 
admissible functions, the main way that it can fail to happen is if the kernel basis functions themselves are 
not linearly independent. These possibilities are not usually addressed in connection with general discussions 
of the minimum norm collocation property, but in many settings linear independence may not always be 
transparent, especially if combinations resulting from linear operators acting on kernel basis functions are 
allowed. (These possibilities can be seen from, among other things, the consequences of the fact that various 
restricted classes of functions, such as polynomials of fixed degree, may have a reproducing kernel.) When 
(|10P is invertible the source parameters (mk) are uniquely determined and in some sense one can say that the 
solution to the inverse point source problem has been obtained. To preserve and extend these inverse source 
interpretational possibilities, the conservative stance is adopted here of requiring that all admitted basis 
function sets be invertible and that the solutions to (fT0|) replicate the specified Ak values. This condition 
is called the generalized least squares collocation (GLLSQC) condition. There are two obvious ways to 
enforce this condition: either on a computational case by case basis or by proving general theoretical results 
about classes of particular basis functions. While it is not comparatively well known outside geophysics, a 
demonstration exists that shows point mass basis functions are independent in R™ for finite n > 1 |23j . Thus, 
in a theoretical sense the GLLSQC condition holds for point mass fits, but on a case-by-case basis some care 
is generally required in solving (|10| due to ill-conditioning. (As a matter of practice, either a singular value 
decomposition or a Hausholder triangulation algorithm implemented to an appropriate number of significant 
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digits should be used.) Finally, it should be readily apparent that while continuous distributions may 
well satisfy the GLLSQC condition, counter-examples can be easily constructed, with perhaps the simplest 
example resulting from the consideration of several concentric homogenous spherical shells. 

With regard to the overall functional analysis setting, the standard course taken nowadays is to adopt 
some type of general Banach space setting (such as one type or another of Sobolev space), where the 
completeness of Cauchy sequences is presupposed. The limit of sequences of functions composed from basis 
functions that satisfy the GLLSQC condition may not satisfy it, hence admitting limits of sequences can 
lead to unwanted interpretational difficulties here. Obviously, to solve (|10[) it is only necessary that a finite 
linear span of basis functions be admitted, which requires only an inner-product structure. In accord with 
the conservative stance outlined above, a structured pre-Hilbert space setting is adopted since a pre-Hilbert 
space setting presupposes an inner product structure, but it does not make the usual assumptions about 
admissible sequences of functions. (This way of specifying functions of interest is somewhat dated, but prior 
to the mid-1950's it was the prevailing way of addressing Dirichlet inner product spaces and it was used, 
for example, in [3]). The adjective "structured" here means that further auxiliary conditions are imposed 
on the class of admissible functions. One such condition is that any set of basis functions considered must 
satisfy the GLLSQC condition stated above. Another requirement is that all functions, /, must be harmonic 
over flj (including dflj). Additional structure is needed to insure that (weighted) Dirichlet integrals over 
unbounded domains can be regarded as defining a positive definite norm. This can be accomplished by 
assuming one (or all) of the following four largely equivalent requirements: 

(I) D[/, /, 1, Slj] is bounded and / tails off at least as fast as l/r as r — ► oo. 

(II) The first Sobolev norm of / over is bounded: D[/, /, 1, Clj] + JJJ \f\ 2 dV < oo. 

(III) / is a potential function generated by a (well-behaved) localized source distribution in Slsj C \ dCl'j. 

(IV) A well behaved series representation for / always exists in the form of a spherical harmonic expansion 
of / in powers of l/r, for the exterior of some sphere contained inside O^-. 

Notice that (I) is the historical approach [3] and that (II) equivalent to (I) for Slj = f^o since both the 
integrals occurring in the first Sobolev norm must be separately bounded. Also for all / not identically 
zero, < D[/, /, /io, Oo] < D[/, /, 1, Clo]/Ro, where positivity can be easily proved using several standard 
properties of harmonic functions. (Specifically, if |/| ^ at some point in ilo then taking the line integral of 
V/ from this point to a point at infinity one can infer that |V/| > for at least one point along this line. 
Then from the mean-value theorem and the maximum-modulus theorem d ■ V/ > must occur throughout 
some neighborhood near dQo for one fixed direction d or another and from this it follows immediately that 
D[/, /, /io, Oo] > 0.) It is clear that (III) and (IV) are largely equivalent to (I) or (II) [although (III) and 
(IV) do not require harmonicity as a separate condition]. Requirement (III) is more or less tantamount to 
the assumption that a collection of a point sources is being modeled. 



3 Half-Space (Qi) Relationships 

For the class of admissible functions just described with bounded source region, the following half-space 
analog of Poisson's solution to the Dirichlet boundary value problem holds p. 268]: 

(12) W(X) = ^f°° f°° m X ",y",0) Z d X »dy» 



^J-oo J-oo [(a: - x") 2 + (y - y") 2 + z 2 } 3 / 2 
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Observation Type 


Source Type 


Invariance Property 


Potential 


Point Mass 


Scalar 


Gravity 


Point Dipole 


Vector 


Gravity Gradient 


Point Quadrupole 


2nd Rank Tensor 



Table 1: Point Data/Source Correspondences for R 3 Half-Space 



1 

8^ 



Vll 1 ■ VW dV 



1 

'8tt 



oo /• oo 



w 



a 



where z > and X" E dfti. Using the facts that the surface integral at infinity is zero for the unbounded 
region f2i, that d / d n := — d /dz and that dS — dxdy for the boundary plane (<9f2i) and applying Green's 
first identitj0 yields: 

(13) {r k \ W) Ei 

Next observe that 
(14) 



— oo J—oo 



dz \X-X'\ 



dx dy 



2 = 



dz \X-X' k \_ 



z=0 i(x-x' k ) 2 + (y-y' k ) 2 + (4) 2 } 3/2 1 

where z' k < 0. Combining (fT2"|). ([T5]) and (T2]) produces: 

(15) W) Ei = W(x' k , y' k , -4)/4 , 

which shows (U}. This relation can be used to evaluate the terms Tk : k' and A k occurring in (fTO"]) : 



(16) 



Tk,k' 



1 



4v / (4-40 2 + K-^) 2 + (4 + 4-') 5 



wk, i4i) 



Two final points are relevant. First, dipole and other higher order multipoles also can be easily fit. In 
general potentials for these point sources can be written as ^2 k Skihkiih -1 ) where Ski are the associated 
source strengths and the Hki are appropriate linear differential operators that can be expressed as a sum of 
partials of various orders with respect to x, y, or z. For example, an electrostatic or point mass dipole term is 
proportional to Dk ■ Vi k . Since Vf^ 1 = — V' k £ k , where V' k :— (d /dx' k , d /dy' k , d /dz' k ) T , and components 
of X' k serve only as parameters when they occur inside the inner product here, all X dependent derivative 
factors operating on l k x that occur inside inner products can be replaced with X' k derivative factors; hence, 
these differential operators can be moved inside or outside the inner products entirely as desired so that 
all required inner products for Tk,k> and Ak can be easily evaluated in closed form. For half-space, these 
possibilities and the associated measurable point quantities are displayed in Table [1] Analogous possibilities 
exist for the spherical exterior case. 

Second, to improve the condition number of the matrix T in (|10p . it is often desirable to use normalized 
basis functions (pk in place of i k : 



(17) 



N k 



1 



0k := — , where Nk '■= 

* 4 ' lie 1 ! 



1 fffnW* + W ■ V<« dV = ff an 4,% dS 
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iii the general norm setting. Introducing T k ^ := (tp k , Cp k ,) = N k N k ,T ktk >, A k := (W, ip k ) = N k A k and 
fh k := m k /N k allows (fTo]) to be reexpressed as 

(18) Tk > k ' ™ fe ' = Ak ■ 

k'=l 

For the R 3 half-space energy norm 

(19) ^ = 2*^1, T ktk ,= . . ==g^^=_ _= ; 4 = ^f?(4i-4), 

v v K - 4' ) + (Vk - y k > ) + (4 + 4 ) v 2 

so that (|18[) can easily be inverted to determine the values of fh k and thus m k . While the use of normalized 
basis functions is not always required for point mass fits, their use in mixed type point source fits is always 
highly recommended due to the diverse associated physical scales that occur there and the attendant large 
condition numbers for T . (An experiment was performed on the results from one of the global NLLSQ 
combined point mass/dipole fits discussed in Section [5l At the specified source locations a linear fit was 
done both with and without normalized basis functions. The ratio of the two resulting condition numbers 
was over 10 20 .) 



4 Spherical Exterior (Qq) Relationships 

The goal of this section is the derivation of the two relevant spherical exterior DIDACKS relationships, ((2^) 
and ((IK). Matters are more complex for this case than they were for the half-space case, but all of the 
supporting issues raised in Sections [2] and [3] can obviously be carried over here, so they are not repeated. 
Let / and g be two admissible functions as discussed in Section [2] and consider the integral norm [20] : 

?2 r r e>2 r r g 

da , where D r := 



(20) (f ) g\~ ~^JJ V r (rfg)da = -^JJ [V r (frg)] 



=R ' dr 



The last expression on the RHS here follows from the evaluation convention of (|4]) . The label "integral norm" 
was chosen by the author since, as discussed below, the integrals required for point mass fitting in f^o can 
be evaluated in closed form and there is little chance of confusing this norm with the usual norm of square 
integrable functions. 

Applying Green's first identity to the f2 energy inner product and noting that the surface integral at 
infinity vanishes, while dS = R^da and d /dn := —T) r for the bounding inner exterior surface of fig, yields: 

(21) (/, ,g) Eo := ^- III V/ • V 5 dV = -ft ffg n rf da = _i| ff [D r {fg)\ 



From (|2U)) and (|2"Tj) it follows immediately that 

(22) (f,g\ = 4R (f,g) Eo - R 2 (f,g) a . 

After the second DIDACKS relation, (J2^), is addressed, it will be shown that the integral norm is proportional 
to the weighted Dirichlet integral, which will complete the proof of the third DIDACKS relationship, (|3Ji). 
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First, observe that the following two equations can be shown through a relatively straightforward evalu- 
ation of their respective LHS and RHSs: 



(23) 

and 
(24) 



r=Ro 



1 



\x-P k \ 



r=R 



Ik 



r=Ra 



D, 



\X- P k I 



r=Ro 



where P k is given by J2b) with X' k = X' and r' k := \X' k \. Employing and yields 



(25) 




d<7 = — 



r=iin 




(2W) 

I^-Pfel 



- WD r 



\X-P k 



Recall that W is harmonic for r > Rq. Applying Green's second identitjH to (f2"5")) yields 




da 



r=R 



Pq 
47rr fe 




\X- Pu 



dV . 



fin 



G?<7 . 



Then using V 2 ( — Pk\ ) = — 47r6(X — Pk), where S is the Dirac delta function [TUl p. 35], gives 



(26) 




da = ^W[P: 



r=R 



or finally, with Pk = \Pk 
(27) 



P k W 



Various other ways to prove (|27| exist. One way is to substitute spherical harmonic expansions for 
W and i k l into the LHS of (|26|) . A second way is to first observe that , ^T/ )j can be evaluated in 
closed form by using spherical coordinates on the LHS of (|26|) . Next, this result can be generalized by 
substituting the integral form of Poisson's equation (not Poisson's integral) for W into the LHS of (f2"7) and 
then reintroducing the closed-form expression just obtained for (t^ 1 , )r > which results in the RHS of ([27)) 
reexpressed in terms of the integral form of Poisson's equation. As previously noted these various steps were 
the ones originally followed by the author and account for the nomenclature. Other proofs have also been 
discovered. 

This leaves the proof of the relationship between the weighted Dirichlet integral and the integral norm. 
Recalling the limit convention implicit for integration over a, (|4]), and expanding the RHS of ([20)l yields 



(28) 



E>2 



rVr(fg) da 



He 



on 



dS = ///n^V^-^V 2 *) 



dv 
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Next temporarily ignore the common factor of R^/Att and consider the two integrals on the RHS of 
Using the fact that D r f = (X/r)-V/, the first integral on the RHS of ([2"5]) can be written as 

dj fgda = J! I I D r {fg)dr Ida = JJJ r- 2 V r (fg) dV 

(29) a a v r =Ro ) 

= 111(9 r~ 3 ) (X-Vf) dV+ffjif r- 3 ) (X-Vg) dV . 

O O 

Green's second identity in the following form will be useful in reexpressing the second term on the RHS 

of 

(30) JJJS7ip-S7cj)dV = -JJr 2 iP(T> r <f))d<7 

Oo f 

where (f> is harmonic in Qq, but ip need not be and both must vanish sufficiently fast as r — > oo. Using this 
identity at the appropriate place twice in the following expression (first with ip = f/r and <f> = g and then 
with / and g reversed), the second term on the RHS of (12"5]) can be rewritten as 

(31) - // rf(V r g)da- ff rg{T) r f)da= fff V<?-V(//r) dV + fff V/-V(<?/r) dV . 



Finally using V(//r) = r _1 V/ - r^Xf in ([31]) and substituting this result along with ([29]) into ([28]) 
produces 

(32) (f, g)^^!!! 'r-'Vf-VgdV . 

Oo 

With the aid of {22]), {27} and ([32]), the second and third DIDACKS relationship [i.e., (0i) and J3ji)] follow 
immediately. From (|32|) and the discussion at the end of Section [2] the integral norm is positive definite. 

DIDACKS dipole and other higher order multipolc implementations differ for the half-space and spherical 
settings in one significant way. While closed-form expressions for all the required inner products for an 
integral-norm point-mass fit can be evaluated just as discussed in Section [3] for higher-order multipole fits 
all the derivatives of the potential for all the lower orders are also required in the spherical case because 
taking partials with respect to the components of XL yields additional terms on the RHS of (|2"T]) . This means, 
for example, that a dipole fit requires not only point gravity information, but point potential information as 
well. It is thus natural in this case to perform not only a point dipole fit, but a combined point mass/dipole 
fit. With this understanding the spherical case analog of Table [T] is readily obtained. 

When an integral norm higher order multipole fit is desired and the lower order derivatives of W are not 
available then it still may be possible to do the fit [20]. For example, if a dipole fit is desired and information 
about W is not available, but the gradient of W is known along various intersecting lines, then potential 
data in the form IU(X) = Wo + 5W(X), where Wo is an unknown constant, can be assumed anywhere 
along these lines since the form of 8W(X) can be found by numerical line integration. For any assumed 
trial value of Wo, a DIDACKS fit can be performed. The results of this fit can then, in turn, be substituted 
into a standard LLSQ type of cost function that is the analog of ([3"3"]) and is based on minimizing gravity 
computations at various sample points. If gravity data itself is plentiful, then an outer-loop optimization 
process can be based on minimizing this new cost function, where Wo is treated as an unknown NLLSQ 
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parameter. In this outer-loop process sample point gravity differences are minimized throughout the fit 
region of interest (which may be only a small part of f2o). When its choice is not obvious, i?o can also be 
treated as a parameter and optimized in the same fashion. 

5 Broader Point Mass Fitting Context 

The goal of this section is to illuminate certain aspects of relevant geophysics and point mass fitting back- 
ground material by way of a brief synopsis (as such this section is somewhat subjective). Before considering 
the associated geophysical context it is useful to clarify the distinction between NLLSQ and GLLSQ or LLSQ 
problems and to briefly consider how the pertinent DIDACKS applications history fits into these categories. 
In the sequel, the distinction between GLLSQ and LLSQ is generally dropped and only the acronym LLSQ is 
used so that both classes can be referred to jointly. If all the source locations are known then linear equation 
sets result for the source strengths and in the DIDACKS approach these linear equation sets are exact due to 

and ([2^) ■ Alternatively, if both the locations and strengths are to be determined then a NLLSQ problem 
results since source locations enter as non-quadratic parameters in the cost function. While only LLSQ prob- 
lems have been addressed so far in this article, accurate low degree and order spherical harmonic (tesseral) 
NLLSQ fits have been obtained by the author and this was, in fact, the first area of DIDACKS applications 
in the early 1980's [20]. These NLLSQ fits are discussed further below. While rather varied approaches 
have been used by different researchers for point mass based geoexploration inverse source applications, for 
the gravity modeling and estimation problems dealt with here far fewer point mass based approaches have 
been employed; nevertheless, not only have the associated research efforts been internationally diverse, but 
the corresponding literature is also extensive so that only a small part of it can be considered here (see [55] 
for additional history and references). LLSQ and NLLSQ point mass gravity models can also be produced 
to serve as synthetic gravity models with realistic attributes, but DIDACKS applications in this area have 
not been considered so it is not discussed below. Finally, two conventions are adopted in the sequel: (a) a 
spherical harmonic expansion to degree and order N will be called a N x N field or expansion and (b) both 
field modeling and estimation problems will be referred to as field reconstruction problems. 

Next, consider the relevant geophysical aspects. Approximately a quarter of a century ago one could 
divide the Earth's gravity field into three parts corresponding to their respective data sources: 

(I) Global spherical harmonic field data derived directly from satellite tracking data and historically con- 
sidered accurate to around the degree and order 8 to 12 range. Here this part of the field is called 'low 
degree and order" and it is taken to be 12 x 12 and below. 

(II) An intermediate field taken here to be the part above 12 x 12 and below 120 x 120, which before the 
evolution of more advanced radar equipped satellites could not be accurately determined. 

(Ill) Regional measurements of geophysical surface quantities known as gravity anomaly and vertical de- 
flections. [These data sets contain part (I) and (II) contributions unless they are factored out.] 

Parts (I) and (II) together will be taken as comprising the global part of the gravity field. As discussed 
below, much higher accuracy is desirable for part (I) than part (II) or part (III) data. Currently very 
accurate spherical harmonic expansions to a much higher degree and order are available on the Internet, 
so the distinction between the above three data sets is not as distinct as it once was, but to understand 
the context of the methods discussed in this section it is useful to keep this historical data partition in 
mind. {Expansions beyond 360 x 360 are common and the recent Gravity Recovery and Climate Experiment 
(GRACE) [91 [19] has already established new global gravity accuracy benchmarks up through 110 x 110.} 

Next consider the origins of part (III) data. Regional data measurements are often made at sea and it is 
in this context that the concepts of gravity anomaly and vertical deflection components are best understood. 
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Suppose that the Earth were composed of a homogenous liquid with the same overall mass and volume, then 
due to rotational effects it would take on an ellipsoidal shape. Geophysicists set up a mathematical model of 
this configuration that they call the reference ellipsoid. The associated field is called the normal gravity field 
and it is the predominant part of the field. If transient effects (like tides and ocean waves) are accounted 
for, then under the influence of gravity the oceans form an equal-potential surface (otherwise water would 
flow from one part to another until an equilibrium was reestablished). If normal gravity is subtracted off 
and if these transient effects are properly taken into account then the measurement taken by a gravimeter 
on a ship is called gravity anomaly, Ac/, and it is a scalar since the measurement is taken along a plumb 
line. The angular displacement of this plumb line from the vertical is called vertical deflection, which can 
be resolved into north-south and east-west components. Due to the equipotential effect just noted, this 
measurement is displaced from the specified reference ellipsoid by a distance that is called geoid height, but 
the measurement is recorded as if it had been made at a point on the reference ellipsoid itself Geoid 
height and scalar potential values are connected by Bruns formula [5] • Geoid height itself can be determined 
by a radar equipped satellite that has a known location. Gravity anomaly is typically measured in units of 
milligal, where 1 milligal = 1 x 10 _5 m/s 2 . A 1 milligal error is considered more-or-less acceptable for regional 
gravity anomaly data processing [9] p. 274] and various underpinning geophysical relationships are derived 
with an inherent approximation consistent with this 1 milligal requirement [U [9] . Much higher accuracy is 
desirable for part (I) than part (II) or part (III) data since gravitational effects are cumulative for most uses 
and low degree and order errors tend not to cancel out. GRACE and other modern fields have error levels 
considerably under a milligal for part (I) data and for this part of the field it is also desirable to have point 
mass models with errors somewhat under a milligal. 

Historically, there have been three primary motivations for performing NLLSQ fits to the low degree and 
order spherical harmonic part of the Earth's gravity field: (a) To find a more efficient computational scheme 
for gravity evaluations, (b) To gain some insight into the distribution of matter in the Earth's interior, (c) 
To conduct goal oriented pure research. Due to the computational ease and speed of low degree and order 
spherical harmonic gravity evaluations that has resulted from computer hardware and software advances, 
(a) has long ceased to be a realistic reason for using point mass fits and this point is totally irrelevant now. 
Other measurement programs and advances have also emerged to address the issues raised by (b), but in fact 
it was never clear that the small number of point masses used in NLLSQ fits could provide truly significant 
mantle or deep core density information, which leaves only (c). NLLSQ point mass fitting presents a very 
challenging problem that can conceivably serve as a sort of test bed for developing techniques to tackle 
other ill-conditioned problematic NLLSQ problems; moreover, aside from these NLLSQ aspects, it is an 
intrinsically interesting potential theory problem that has associated cross-fertilization possibilities. 

The first step in attacking any LLSQ or NLLSQ problem is to set up a cost function. A commonly 
chosen minimization philosophy for the NLLSQ point mass fitting problem is that of matching the observed 
quantities as closely as possible. Since there is a classical result in geophysics (Stokes' integral) that says 
that if gravity anomaly is known over the entire reference ellipsoid then the external field quantities can be 
reconstructed, standard approaches to performing global NLLSQ low degree and order point mass fits often 
have been based on matching gravity anomaly at Ni different sample points, Xi, specified on the reference 
surface. This technique will be called the "classical point mass" approach. If AgpM (Xi) denotes the anomaly 
generated by a collection of Nk point masses and Ag re f(Xi) represents the truth anomaly value at the same 
point, this classical NLLSQ point mass approach can be framed through the requirement that the following 
cost function be minimized: 

Ni 

(33) * = £ ( ^9PM(Xi) - Ag ref (Xi)) 2 , 

i=l 

where Ni >> Nk- The usual prescribed number of point masses, Nk, is approximately 50 for a NLLSQ 9x9 
fit, while Nk « 80 for a 12 x 12 fit. The resulting gravity anomaly error standard deviations (sigma) achieved 
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by applying this classical point mass fitting approach has typically been from 3/4 to several milligals. These 
approaches normally have a much more sizeable error in the smaller degree and order terms than is desirable. 

The global NLLSQ point mass fitting problem itself is inherently nonlinear and the point masses must 
be quite deep to obtain good results, which heightens numerical difficulties and associated convergence 
problems. Given this state of affairs, coupled with the facts that the smaller degree and order errors cannot 
be easily removed using this classical approach and that it contains inherent sampling and discretization 
error, there is an error floor of around 2/3 milligal that these approaches historically have not been able to 
overcome. While the associated regional LLSQ approaches that have been tried are remarkably diverse, far 
fewer low degree and order NLLSQ approaches have been employed — nevertheless the diversity of attempted 
NLLSQ point mass approaches in the literature is much wider than the above discussion indicates, but still 
the accuracy levels and deficiencies of the classical point mass approach presented above are thought to be 
representative of these other existing NLLSQ attempts as well. In Section [2] it was argued that the usual 
philosophy of matching measured quantities is not the correct one and that an energy basis or weighted 
energy basis is clearly called for in approaching both LLSQ and NLLSQ gravity modeling problems. 

As noted above, accurate global NLLSQ low-degree and order spherical-harmonic point-mass fits using 
the DIDACKS approach were first obtained almost a quarter of a century ago by the author and in the past 
as better spherical harmonic data has became available more accurate fits have been obtained [20]. This 
has been an ongoing effort and the current NLLSQ 50 point mass fit to the 9x9 part of a recent field has 
a sigma error of about 0.035 milligal, while the corresponding error in an NLLSQ 80 point mass fit to the 
12 x 12 part of the field is 0.030 milligal. As one might expect, additional masses here can offer marked 
improvements in accuracy, but at a loss of efficiency. Combined point mass/dipole fits (c.f., Section [4]), 
which are combinations of point masses and dipoles at the same location, have also been performed over the 
years 20 . Only about half as many of these combined point sources are required for an accurate NLLSQ 
fit: 22 for a 9 x 9 field and 35 for a 12 x 12 field. These NLLSQ DIDACKS fits have all been based on the 
integral norm introduced in Section 2J Besides the basic DIDACKS formalism, these fits also use additional 
specialized NLLSQ techniques developed by the author to handle the existence of numerous false minima at 
various physical scales. 

While the fits should be doable, neither LLSQ nor NLLSQ DIDACKS fits to the complete global inter- 
mediate part of the field [part (II)] have been attempted, although efficient and accurate DIDACKS fits for 
various regional and local gravity modeling and estimation problems using both LLSQ and partial NLLSQ 
implementations have been obtained by the author [20] . Regional point mass gravity field reconstruction 
[part (III)] based on the "classical point mass" approach epitomized by (|33[) has also often been successfully 
done over the years by various researchers; however, some have reported disappointing results. This is not 
surprising since considerable patience is also often required in order to obtain the best or even good results 
with the DIDACKS approach due to source placement issues. Thus while both DIDACKS and the classical 
point mass fitting approach can be viewed as an alternative to GC for regional applications, they both clearly 
share the same sensitivity to and dependence on point mass placement. This sensitivity to point mass place- 
ment is also apparent from the low degree and order NLLSQ DIDACKS results since NLLSQ iterations have 
obviously reduced the errors by several orders of magnitude over that of initial trial configurations. GC is 
designed to not be as sensitive to the placement of the kernel measurement points (which is the corresponding 
placement issue for it) and thus GC generally requires less patience and skill. Alternative part (II) and (III) 
grid-based point mass approaches are referenced and discussed in [25] . 

It is thus useful to briefly describe GC in order to compare and contrast it with the DIDACKS approach. 
Unlike DIDACKS, GC theory generally assumes all available data is used. GC as commonly practiced differs 
from standard collocation techniques utilized by applied mathematicians in five basic ways: 

1 A SRK basis is assumed and emphasis is primarily focused on a statistical (covariance) interpretation 
given to these kernels. 
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2 Laplace's equation in R 3 is always assumed to hold over the region of interest. 

3 This field region is assumed to be either half-space or the exterior of a sphere (but the half-space setting 
is used only occasionally for localized regional distributions, so it is not stressed). 

4 The study of kernels is broken down into the standard empirically modeled statistical covariance kernels 
and analytical collocation (where closed-form kernels are studied) — analytical collocation is generally 
only used when the functional form of the kernel is understood to be a workable representation of the 
statistical covariance function. 

5 Field measurement errors are allowed. 

Practitioners of GC have cataloged most, if not all, useful kernel possibilities allowed by these five differences. 
The main interpretational basis of GC is the GC or minimum norm property previously mentioned in 
connection with the GLLSQC condition in Section [Jl Under the assumption of errorless measurements, 
from all possible candidate functions that reproduce the given point measurements, GC selects the one that 
is smoothest — that is the one with the smallest norm, where the norm is determined by the underlying 
covariance function [18[ pp. 207-220]. It this minimum norm property, in concert with the emphasis on 
statistical covariance functions indicated by the first point above, that makes GC approaches less sensitive 
to kernel point placement issues; however, as indicated above this also implies a corresponding loss of fitting 
responsiveness and thus, for example, it would be hard to argue that GC is capable of the overall level of 
economy and efficency indicated above for the NLLSQ DIDACKS low degree and order tesseral fits. Further 
observe that just as in RKHS theory, GC also satisfies a least squares norm property [18], but in practice 
this property is usually ignored by physical geodesists. 

For LLSQ gravity reconstruction problems GC has historically been the first procedure of choice and 
for these problems one could argue that when good covariance data is available, these techniques are safe 
and easy to use; however, as classically practiced GC techniques do have notable limitations, which are not 
shared by the DIDACKS approach: 

(A) NLLSQ applications cannot be treated. 

(B) GC techniques are customarily applicable only to the Earth's gravitational field environment since they 
require covariance information that is often only available in this context. 

(C) The application of GC requires a certain level of familiarity and thus it is almost never adapted 
to problems outside the geophysical realm, even when appropriate covariance data can be gathered. 
(There are, however, various other areas that use somewhat related kernel techniques [21].) 

(D) Inverse source estimation problems cannot be entertained. 

In summary, for gravity reconstruction problems with accurate selected data sets where either approach 
can be used, results using the DIDACKS approach can be both much better or much worse than one 
might normally expect with GC approaches since the ingenuity, implementation skill and patience of the 
practitioner have a much more significant bearing on the outcome for DIDACKS approaches. GC is and 
probably will always remain the primary mathematical tool employed for raw gravity data processing and 
related uses where efficiency is not the primary concern, since its behavior in these arenas is well understood 
and it can handle measurement errors naturally. 

6 Mathematical Connections to Geophysical Collocation 

As indicated in Section [TJ ([2^,) can be related to a line of preexisting research performed by Krarup in 
conjunction with his studies of GC [14] . Also an independent parallel line of point source research exists 
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that maintains direct connections to GC itself and can be considered a direct outgrowth of Krarup's original 
work. This alternative collocation based point source scheme is briefly considered after the connections of 
DIDACKS theory to Krarup's work are addressed. Since this alternative scheme and Krarup's work are 
primarily based on GC for spherical exteriors, only this geometry will be considered in the sequel. 

Krarup first introduced the weighted integral occurring on the RHS of f|32|) in conjuction with his study 
of GC [TJJ pp. 62-65]. The SRK corresponding to this norm (the "Krarup kernel," K K ) has the form 

(34) K K (P,Q)-- H " 



RI-2P-Q + (PQ/R ) 2 



which Krarup and his followers extensively studied and applied. In (|34| Q = \Q\- To understand how point 
mass fitting enters, observe that this SRK can be recast in the form R 2 /(P\Q — X'\), or |X'|/|Q — X'\ 
since X' = PRq/P 2 , which is proportional to a point mass potential at a fixed location. Making the ansatz 
X 1 — > X' k and Q — > X allows one to transform a collocation fit into a point mass fit if the mass locations 
are restricted to be at a fixed depth r' k = \X' k \ — constant, where r' k is simply treated as an overall constant of 
proportionality. The relevant history of this and related ideas is briefly touched on below. Here one possible 
next step, which Krarup and his followers did not apparently take since it entails abandoning symmetric 
kernel forms altogether, is to generalize this procedure to independent variable depths by absorbing the factors 
Pk '■= Ro/r' k and Rq into each of the collocation fitting parameters separately and then reinterpreting the 
resulting collocation fit as a point mass fit, which yields a fit based on (J2ji) as the end product. 

Conversely, DIDACKS point mass fits based on the integral norm can be reinterpreted as collocation 
fits [20]. The resulting fits are based on what is called the reciprocal distance covariance function [THJ p. 
182], which corresponds to Rq/PQ times the Krarup kernel specified by (|34[) . Notice that while the integral 
norm can be reexpressed in terms of the "Krarup norm" by (|32|) , it is useful to retain the integral norm as 
a distinct entity specified by ([20]) since (a) it is a surface integral rather than a volume integral, (b) there 
are some consequences of this form, such as (|22|) , that are not apparent from the weighted Dirichlet integral 
form itself and (c) the Krarup norm is associated exclusively with the Krarup kernel in SRK form and it is 
primarily linked to GC theory where, as previously noted, the goals and practices are quite different. 

As summarized in T5] and as just indicated, starting in the early 1980's a distinct line of point source 
research based on Krarup's original work above was developed by A. N. Marchenko and others [7J. This 
research is based on maintaining connections to symmetric global GC covariance kernel forms. For the 
spherical case, assuming the usual GC covariance properties, the general form of allowed global covariance 
kernel, C(P, Q), can be written as: 

(35) C(P, Q) = V k n — -M P„(cos^), where cos?/> = (P-Q/PQ), 



=o 



where the P„ are standard (un-normalized) Legendre polynomials and the k n are constants [18| p. 181]. 
Here, as before, i?o is the radius of the spherical region. Mathematically, a kernel specified by ([33)) can, 
in general, be simply considered a SRK with an added layer of statistical interpretation. Solving <[2t>) for 
X 1 = X'(P) and substituting the result into ([35]) yields a kernel that is a function of one interior point and 
one exterior point. For certain applications, when (|35|) is reexpressed in terms of a smaller radius Rb < Ra 
the resulting sphere is known as a Bjerhammar sphere. For particular choices of k n if (|35p can be rewritten 
as a closed- form expression and if this expression has the correct form, such as (|34[) . then the result can be 
reexpressed as a point mass or other linear combinations of point or line source potentials. As just noted 
in connection with ([34]) . additional parameters will occur in the resulting expressions that do not appear in 
the source potentials themselves, but if the locations of the sources are restricted these parameters can be 
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explained away as common constant factors. A Bjerhammar sphere is generally used since it allows for an 
independent adjustment of the overall depth of the interior source points. By maintaining connections to 
the symmetric kernel forms that are allowed by (|35[) . a statistical covariance interpretation is possible. 

As one might expect there are notable differences between this approach and DIDACKS since this al- 
ternative approach adopts the basic philosophy of maintaining connections to GC and DIDACKS does not. 
In particular for a specified domain, from RKHS theory one can deduce that the choice of inner product 
(or norm) and reproducing kernel must be in one-to-one correspondence [18j ; hence, the choice of multipole 
form to be fit and norm are directly linked. A list of known covariance kernel and multipole point source 
correspondences for this approach can be found in [15) . For example, in [7J some of the consequences of 
the Krarup norm choice, along with its half-space approximation, are examined for point mass fits. By 
employing (|35p this approach allows for an underpinning statistical covariance interpretation, which is im- 
portant for many geophysical applications; nevertheless, implicitly retaining symmetric kernel forms also 
adds a layer of additional complexity, which has the effect of greatly complicating the theory and of limiting 
the possible choices of dipole and higher multipole forms. For uses outside geophysics there are also clearly 
interpretational difficulties that limit its use. 

In DIDACKS theory all attempts at retaining connections to SRKs are dropped and the primary empha- 
sis is placed on the (weighted) energy or integral norm and the associated kernel form As noted earlier, 
this makes consideration of dipoles and higher multipoles trivial since the required inner products are readily 
obtained. Moreover, since the DIDACKS approach maintains the same norm choice for all source types there 
are few interpretational issues — especially with regards to multipole fits of all orders. The types of geophys- 
ical areas where one of these two approaches or GC should be preferred over the others clearly warrants 
further study and consideration since, to date, there have been no researchers proficient in applying all three 
algorithms. (While GC has had many practitioners and this alternative collocation based approach has had 
a few practitioners, so far DIDACKS development and use has been limited to the author's involvement.) 
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